**********************************************************************
c                         SWIFT_SYMBA7_HYBRIDGR.F
c**********************************************************************
c
c This is a version deviced for study of terrestrial
c in the 5-planet scenario. Propagates only planets are read
c from an input file provided by Symba integration. Planet (nplalost)
c in the sequence may be lost and the code takes care of it.

      include './swift.inc'
      include './hybrid.inc'

      real*8 mass(NTPMAX),j2rp2,j4rp4
      real*8 xh(NTPMAX),yh(NTPMAX),zh(NTPMAX)
      real*8 vxh(NTPMAX),vyh(NTPMAX),vzh(NTPMAX)
      real*8 axh(NTPMAX),ayh(NTPMAX),azh(NTPMAX)
      real*8 vxb(NTPMAX),vyb(NTPMAX),vzb(NTPMAX)

      real*8 masw(NTPMAX),xhw(NTPMAX),yhw(NTPMAX),zhw(NTPMAX)
      real*8 vxhw(NTPMAX),vyhw(NTPMAX),vzhw(NTPMAX)
      real*8 rpw(NTPMAX),rhw(NTPMAX)

      integer i1stv,i1sta,ierr
      integer nbod,ipl,ialpha,nstp,nbodo,ntot,nbodm
      integer iflgchk,iub,iuj,iud,iue,ium,id
      real*8 dtp,a,e,inc,capom,omega,capm

      real*8 t0,tstop,dt,dtout,dtdump
      real*8 t,tout,tdump,tfrac,tpl,tmpl,dr

      real*8 rmin,rmax,rmaxu,qmin
      real*8 rpl(NTPMAX),rhill(NTPMAX),mtiny,eoff
      integer isenc,mergelst(2,NTPMAX),mergecnt
      integer iecnt(NTPMAX)
      logical*2 lclose

      character*80 outfile,fopenstat,inplfile
      character*80 inparfile,orbitfile

      integer io_read_hdr_r,io_read_line_r,io_read_mass_r
      integer i,j,nplalost,ibar

      real*8 xbeg(NPLMAX),ybeg(NPLMAX),zbeg(NPLMAX)
      real*8 vxbeg(NPLMAX),vybeg(NPLMAX),vzbeg(NPLMAX)
      real*8 xend(NPLMAX),yend(NPLMAX),zend(NPLMAX)
      real*8 vxend(NPLMAX),vyend(NPLMAX),vzend(NPLMAX)
      real*8 gm,clight,vbreak

      integer njov,njovo,njovm,nlefj,npl0
      real*8 rpjov(NPLMAX),dth,mjov(NPLMAX),mjbeg(NPLMAX)
      real*8 xtmpj(NPLMAX,0:NINTERP),ytmpj(NPLMAX,0:NINTERP)
      real*8 vxtmpj(NPLMAX,0:NINTERP),vytmpj(NPLMAX,0:NINTERP)
      real*8 ztmpj(NPLMAX,0:NINTERP),vztmpj(NPLMAX,0:NINTERP)
      real*8 mtmpj(NPLMAX),vxbj(NPLMAX),vybj(NPLMAX),vzbj(NPLMAX)

      integer nmaxgr
      real*8 gmc2
      integer tini,tfin

      integer ijob
      character*80 filediscard,filedumppl
      character*80 filedumpparam,massoutfile     
      common /discard/ filediscard

c...  Common block to be passed to symba7_merge_split
      real*8 vbreak2
      common /BREAK/ vbreak2

c...  OMP stuff
!$    logical OMP_get_dynamic
!$    integer nthreads,OMP_get_max_threads

c... Executable code

c...  OMP stuff
!$    write(*,*) 'Dynamic thread allocation: ',OMP_get_dynamic()
!$    nthreads = OMP_get_max_threads()             ! In the *parallel* case
!$    write(*,'(a)')      ' OpenMP parameters:'
!$    write(*,'(a)')      ' ------------------'
!$    write(*,'(a,i3,/)') ' Number of threads  = ', nthreads 

c... get initial time
c      tini = time8()
c      write(*,*) ctime(tini)

      dr = 180.d0/PI

c... Get data for the run and the test particles
      write(*,*) 'Enter name of parameter data file : '
      read(*,999) inparfile
999   format(a)
      call io_init_param(inparfile,t0,tstop,dt,dtout,dtdump,
     &         iflgchk,rmin,rmax,rmaxu,qmin,lclose,outfile,fopenstat)

c... forces output and dump to be at least 1 year (according to ephemeris file)
      if (dtout.lt.1.d0) dtout = 1.d0
      if (dtdump.lt.1.d0) dtdump = 1.d0
      write(*,*) 'New dtout and dtdump',dtout,dtdump     

c... Prompt and read name of planet data file
c... ONLY TERRESTRIALS AND PLANETESIMALS - TERRESTRIALS FIRST !!!
      write(*,*) ' '
      write(*,*) 'Enter name of planet data file : '
      read(*,999) inplfile
      call io_init_pl_symba(inplfile,lclose,iflgchk,nbod,mass,
     &     xh,yh,zh,vxh,vyh,vzh,rpl,rhill,j2rp2,j4rp4)

c... Read the mass threshold for planetesimals mass and breakup veolocity
c... Set vbreak >> 1 (e.g. 10^8) for no splitting
c... Velocity units must be compatible with IC units
      read(*,*) mtiny,vbreak
      vbreak = 1.0d8
      vbreak2 = vbreak**2

c Specify index of the job (<999)
        read(*,*)ijob

 78     format('discard.out_',i5.5)
        write(filediscard,78)ijob
 80     format('dump_pl.dat_',i5.5)
        write(filedumppl,80)ijob
 82     format('dump_param.dat_',i5.5)
        write(filedumpparam,82)ijob
 83     format('output.bin_',i5.5)
        write(outfile,83)ijob
 84     format('mass.output.bin_',i5.5)
        write(massoutfile,84)ijob

c... Calculate the location of the last massive planet (only terrestrials)
      call symba7_nbodm(nbod,mass,mtiny,nbodm)
      write(*,*) 'Of',nbod,' terrestrials',nbodm,' are massive'

c... Prompt and read name of planet data file
      write(*,*) 'Enter planetary output data : '
      read(*,999) orbitfile
      call io_open(90,orbitfile,'old','UNFORMATTED',ierr)
      if (ierr.ne.0) then
         write(*,*) 'Problem in opening planet.bin'
         call util_exit(1)
      endif
      call io_open(91,'mass.'//orbitfile,'old','UNFORMATTED',ierr)
      if (ierr.ne.0) then
         write(*,*) 'Problem in opening mass.planet.bin'
         call util_exit(1)
      endif

c... Specify which planet got lost during the evolution (0 if none)
      write(*,*) 'Enter planet that is lost (from 1 to 5, only jovians 
     & without Sun) : '
      read(*,*) nplalost
      if ((nplalost.lt.0).or.(nplalost.gt.5)) then
         write(*,*) 'Something fishy in specifying the lost planet; ',
     &   nplalost
         call util_exit(1)
      endif
c      npl0 = nplalost+1   ! select this if want to turn off mass of fifth pl
      npl0 = -1

c... Reads speed of light and max planets (from inner to outer) to be corrected by relativity
c...   clight must be in units consistent with those used in initial conditions
c...   set nmaxgr = 1 for no relativistic corrections
      read(*,*) nmaxgr,clight
      if (nmaxgr.ne.1) nmaxgr = nbodm        ! apply GR to all massive bodies !!!!
      gmc2 = mass(1)/clight/clight           ! determines GM/c^2

c... Initialize initial time and times for first output and first dump
      iub = 20
      iuj = 30
      iud = 40
      iue = 60
      ium = 21

c... Radii of planets: note R_Uranus ~ R_Neptune ~ R_5th_planet
      rpjov(2) = (4.772e-4)
      rpjov(3) = (4.011e-4)
      rpjov(4) = (1.698e-4)
      rpjov(5) = rpjov(4)
      rpjov(6) = rpjov(4)
        
c... Reading the first planet data point from planet files
c... NOTE: the mass of the Sun is also read from mass.planet file
 77   ierr = io_read_hdr_r(90,tpl,njov,nlefj) 
      if (ierr.ne.0) goto 2
      ierr = io_read_mass_r(tmpl,njovm,mjov,91)
      if (ierr.ne.0) goto 2
      if ((tpl.ne.tmpl).or.(njov.ne.njovm)) goto 2

      if (nbod.eq.1) mass(1) = mjov(1)

c... reading elements and converting to x,v 
      do i=1,njov
         mjbeg(i) = mjov(i)   ! stores the masses at the begin of the step to be
                              ! used for interpolation later
      enddo
      do ipl=2,njov
         ierr = io_read_line_r(90,id,a,e,inc,capom,omega,capm) 
         if(ierr.ne.0) goto 2
         
         ialpha = -1
         gm = mjov(1)+mjov(ipl)
         call orbel_el2xv(gm,ialpha,a,e,inc,capom,omega,capm,
     &              xbeg(ipl),ybeg(ipl),zbeg(ipl),
     &              vxbeg(ipl),vybeg(ipl),vzbeg(ipl))
c         call flush()
      enddo

c... check if times are consistent
      if (tpl.lt.t0) goto 77    ! necessary when simulation is restarted
      njovo = njov
      t = tpl
      t0 = tpl
      if (tstop.le.t0) then
         write(*,*)' Fishy sequence of t0, tpl, tstop; ',tpl,tstop
         call util_exit(1)
      endif
      tout = tpl + dtout
      tdump = tpl + dtdump

c... Do the initial io write, stores the data such that:
c...  Jovian planets (from interpolation) come first
c...  Terrestrial planets come second
c...  Terrestrial planetesimals come third
      ntot = nbod+njov-1
      masw(1) = mass(1)
      xhw(1) = xh(1)
      yhw(1) = yh(1)
      zhw(1) = zh(1)
      vxhw(1) = vxh(1)
      vyhw(1) = vyh(1)
      vzhw(1) = vzh(1)
      rpw(1) = rpl(1)
      rhw(1) = rhill(1)
      do i=2,njov
         j = i
         masw(j) = mjov(i)
         xhw(j) = xbeg(i)
         yhw(j) = ybeg(i)
         zhw(j) = zbeg(i)
         vxhw(j) = vxbeg(i)
         vyhw(j) = vybeg(i)
         vzhw(j) = vzbeg(i)
         rpw(j) = rpjov(i)
         call util_hills1(mjov(1),masw(j),xhw(j),yhw(j),zhw(j),
     &        vxhw(j),vyhw(j),vzhw(j),rhw(j)) 
      enddo
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) PRIVATE(i,j)
!$OMP&   SHARED(nbod,njov,mass,xh,yh,zh,vxh,vyh,vzh,rpl,rhill,
!$OMP&          masw,xhw,yhw,zhw,vxhw,vyhw,vzhw,rpw,rhw)
      do i=2,nbod
         j = i+njov-1
         masw(j) = mass(i)
         xhw(j) = xh(i)
         yhw(j) = yh(i)
         zhw(j) = zh(i)
         vxhw(j) = vxh(i)
         vyhw(j) = vyh(i)
         vzhw(j) = vzh(i)
         rpw(j) = rpl(i)
         rhw(j) = rhill(i)
      enddo
!$OMP END PARALLEL DO

c... only to grab initial conditions of all the system
c      write(12,*) ntot
c      write(12,*) masw(1)
c      do j=2,ntot
c         write(12,*) masw(j),rhw(j),rpw(j)
c         write(12,*) xhw(j),yhw(j),zhw(j)
c         write(12,*) vxhw(j),vyhw(j),vzhw(j)
c      enddo
c      stop

c... ATTENTION: routine modified here to compute terrestrial elements
c... using mass(1) and jovian elements using mjov(1)
ccc      if (btest(iflgchk,1))  then ! bit 1 is set
         call io_write_frame_r_symba(t0,ntot,nbodm,njov,masw,
     &           mjov(1),xhw,yhw,zhw,vxhw,vyhw,vzhw,outfile,
     &           iub,fopenstat)
         call io_write_mass_r(t0,ntot,masw,outfile,ium,fopenstat)
ccc      endif

c... check for total energy
      eoff = 0.0d0
      if(btest(iflgchk,2))  then ! bit 2 is set
         call anal_energy_write(t0,ntot,masw,j2rp2,j4rp4,xhw,yhw,zhw,
     &        vxhw,vyhw,vzhw,iue,fopenstat,eoff)
      endif

c... must initize discard io routine
ccc      if (btest(iflgchk,4))  then ! bit 4 is set
         call io_discard_mass(0,t,0,mass(1),rpl(1),xh(1),yh(1),zh(1),
     &        vxh(1),vyh(1),vzh(1),iud,-1,fopenstat)
ccc      endif

c*************** here is the big loop *************************************
      write(*,*) ' ************** MAIN LOOP ****************** '
c      call flush()

      do while ( (t.lt.tstop) )

c... Read the next data point in planet files and interpolate
c... first read the mass.planet file
c... NOTE: the mass of the Sun is also read from mass.planet file
         ierr = io_read_hdr_r(90,tpl,njov,nlefj) 
         if (ierr.ne.0) goto 2
         ierr = io_read_mass_r(tmpl,njovm,mjov,91)
         if (ierr.ne.0) goto 2
         if ((tpl.ne.tmpl).or.(njov.ne.njovm)) goto 2
 
         if (nbod.eq.1) mass(1) = mjov(1)

c... check if planet has been lost; WATCH: we said it should be (nplalost)-th planet
         if (njov.lt.njovo) then
            do ipl=(nplalost+1),njov            ! shift indexes of Uranus and Neptune
               mjbeg(ipl) = mjbeg(ipl+1) 
               xbeg(ipl) = xbeg(ipl+1)
               ybeg(ipl) = ybeg(ipl+1)
               zbeg(ipl) = zbeg(ipl+1)
               vxbeg(ipl) = vxbeg(ipl+1)
               vybeg(ipl) = vybeg(ipl+1)
               vzbeg(ipl) = vzbeg(ipl+1)
               rpjov(ipl) = rpjov(ipl+1)
            enddo
            npl0 = -1
         endif
         njovo = njov

c... then read the jovian elements from planet file and convert to x,v
         do ipl=2,njov
            ierr = io_read_line_r(90,id,a,e,inc,capom,omega,capm) 
            if (ierr.ne.0) goto 2

c... the lost planet may be temporarily on hyperbolic orbit!
            if (e.ge.1.d0) then
               ialpha=1
            else
               ialpha=-1
            endif
            gm = mjov(1)+mjov(ipl)
            call orbel_el2xv(gm,ialpha,a,e,inc,capom,omega,capm,
     &          xend(ipl),yend(ipl),zend(ipl),
     &          vxend(ipl),vyend(ipl),vzend(ipl))
c            call flush()
         enddo

c... define the step in planetary ephemeris and substeps
         dtp = tpl-t
         nstp = int(dtp/dt)
c         write(*,*) t,tpl,dtp,nstp,dt
         if (nstp.eq.0)then
            write(*,*) 'nstp equal to zero'
            write(*,*) 'tpl,t=',tpl,t
            call util_exit(1)
         endif
         if (nstp.ne.NINTERP)then
            write(*,*) 'ATTENTION: nstp',nstp,
     &                 ' is diferent from NINTERP',NINTERP
            write(*,*) 'Fix NINTERP in hybrid.inc'
            call util_exit(1)
         endif
         dt = dtp/float(nstp)    ! modify dt so there are EXACTLY nstp in dtp
         dth = 0.5d0*dt

c... Now interpolate planetary positions in dt-steps
c... NOTE: FORWARD propagation for interpolation is done with the planet mass
c...       at the BEGIN of the interval (mjbeg). 
c...       BACKWARD propagation uses masses at the END of the interval (mjov).  
         call rmvs3_interp_mod(njov,xbeg,ybeg,zbeg,vxbeg,vybeg,
     &     vzbeg,xend,yend,zend,vxend,vyend,vzend,dtp,mjbeg,mjov,nstp,
     &     xtmpj,ytmpj,ztmpj,vxtmpj,vytmpj,vztmpj)

c... set temporary masses of the jovian planets (necessary if some
c... mass was set to zero)
         do i=1,njov
            mtmpj(i) = mjov(i)
            if (i.eq.npl0) mtmpj(i) = 0.d0
         enddo

c... Integrate the terrestrial planets
         i1stv = 0
         i1sta = 0
         do i=1,nstp

c... just for safe to avoid passing trash to getacc and lindrift routines
            xtmpj(1,i-1) = 0.0d0
            ytmpj(1,i-1) = 0.0d0
            ztmpj(1,i-1) = 0.0d0
            vxtmpj(1,i-1) = 0.0d0
            vytmpj(1,i-1) = 0.0d0
            vztmpj(1,i-1) = 0.0d0
            xtmpj(1,i) = 0.0d0
            ytmpj(1,i) = 0.0d0
            ztmpj(1,i) = 0.0d0
            vxtmpj(1,i) = 0.0d0
            vytmpj(1,i) = 0.0d0
            vztmpj(1,i) = 0.0d0

c... apply a kick on innermost planets due to relativistic acceleration
c... kick is applied on helio velocities
            if (nmaxgr.gt.1) then
               call getacch_gr(nbod,mass,xh,yh,zh,vxh,vyh,vzh,
     &           axh,ayh,azh,nmaxgr,gmc2)
               call kickvh(nbod,vxh,vyh,vzh,axh,ayh,azh,dth)
               i1stv = 0     ! FORCE TO RECOMPUTE BARYCENTRIC VELOCITIES AFTER GR KICK
            endif

c... if this is the begining of the sequence
            if (i1stv.eq.0) then

c... convert helio velocities to bary, wrt barycenter of the whole system
c... ibar = 0 => input for terrestrial is helio
c... ibar = 1 => input for terrestrial is bary
c... input for jovians is always helio
              ibar = 0
              call coord_vh2b_jov(ibar,nbod,mass,vxh,vyh,vzh,njov,mtmpj,
     &             vxtmpj(1,i-1),vytmpj(1,i-1),vztmpj(1,i-1),
     &             vxbj,vybj,vzbj)
              call coord_vh2b_ter(nbod,vxh,vyh,vzh,njov,
     &             vxbj,vybj,vzbj,vxb,vyb,vzb)

              i1stv = 1    ! NOTE THIS
            endif

            if (i1sta.eq.0) then

c... compute accelerations on terrestrials at the begining of the step
              call helio_getacch_jov(njov,mtmpj,xtmpj(1,i-1),
     &           ytmpj(1,i-1),ztmpj(1,i-1),nbod,xh,yh,zh,axh,ayh,azh)

              i1sta = 1    ! NOTE THIS
            endif

c... apply a kick on terrestrial from jovians for half dt
c... kick is applied on bary velocities
c... use the accelerations calculated at the end of previous step
            call kickvh(nbod,vxb,vyb,vzb,axh,ayh,azh,dth)

c... do a linear drift of terrestrials due to momentum of the Jov for half dt
            call helio_lindrift_jov(nbod,mass,njov,mtmpj,
     &           vxbj,vybj,vzbj,dth,xh,yh,zh)

c... THIS BLOCK CONSIDERS ONLY THE TERRESTRIALS, ALL PERTURBATIONS
c... JOVIANS HAVE BEEN ALREADY ADDED

c... do a normal full step for terrestrials
c... input are helio coordinates and helio and bary velocities
c... close encounters among terrestrial are resolved, but perturbation from
c...  jovians is not taken into account
            call symba7_step_pl_bary(t,nbod,nbodm,mass,j2rp2,
     &           j4rp4,xh,yh,zh,vxh,vyh,vzh,vxb,vyb,vzb,dt,lclose,rpl,
     &           isenc,mergelst,mergecnt,iecnt,eoff,rhill,mtiny)

c... END OF BLOCK WITH ONLY TERRESTRIALS

c... at this point it is necessary to have the barycentric velocities of all bodies
c... terrestrials are already given 
c... jovians must be converted from interpolated heliocentric
            ibar = 1
            call coord_vh2b_jov(ibar,nbod,mass,vxb,vyb,vzb,njov,mtmpj,
     &           vxtmpj(1,i),vytmpj(1,i),vztmpj(1,i),
     &           vxbj,vybj,vzbj)

c... do a linear drift of terrestrials due to momentum of the Jov for half dt
            call helio_lindrift_jov(nbod,mass,njov,mtmpj,
     &           vxbj,vybj,vzbj,dth,xh,yh,zh)

c... compute accelerations on terrestrials at the end of the step
            call helio_getacch_jov(njov,mtmpj,xtmpj(1,i),
     &           ytmpj(1,i),ztmpj(1,i),nbod,xh,yh,zh,axh,ayh,azh)

c... apply a kick on terrestrial from jovians for half dt
c... kick is applied on bary velocities
            call kickvh(nbod,vxb,vyb,vzb,axh,ayh,azh,dth)

c... update helio velocities at the end of the step
c... this is needed for close encounters routines and GR
            call coord_vb2h_ter(nbod,mass,vxb,vyb,vzb,njov,mtmpj,
     &             vxbj,vybj,vzbj,vxh,vyh,vzh)

c... apply a kick on innermost planets due to relativistic acceleration
c... kick is applied on helio velocities
            if (nmaxgr.gt.1) then
               call getacch_gr(nbod,mass,xh,yh,zh,vxh,vyh,vzh,
     &           axh,ayh,azh,nmaxgr,gmc2)
               call kickvh(nbod,vxh,vyh,vzh,axh,ayh,azh,dth)
            endif

c... update the time
            t = t + dt

ccc            if(btest(iflgchk,4))  then ! bit 4 is set

c... check for discards or merges among terrestrials
               nbodo = nbod
               call discard_massive_hybrid(t,dt,nbod,mass,xh,yh,zh,
     &              vxh,vyh,vzh,njov,mtmpj,xtmpj(1,i),ytmpj(1,i),
     &              ztmpj(1,i),vxtmpj(1,i),vytmpj(1,i),vztmpj(1,i),
     &              rmin,rmax,rmaxu,qmin,lclose,rpl,
     &              rhill,isenc,mergelst,mergecnt,iecnt,eoff)

c... calculate the location of the last massive planet (only terrestrials)
               if(nbodo.ne.nbod) then
                  call symba7_nbodm(nbod,mass,mtiny,nbodm)
c                  call flush()
                  if (nmaxgr.ne.1) nmaxgr = nbodm  ! for GR
                  i1stv = 0
                  i1sta = 0
               endif

ccc            endif

         enddo ! here ends the cycle over i

c... Update time according to jovian planets file (before performing check for output)
         t = tpl

c... Update array for writing and dump
c...  Terrestrial planets come first
c...  Jovian planets (from interpolation) come second
c...  Terrestrial planetesimals come third
         ntot = nbod+njov-1
         masw(1) = mass(1)
         xhw(1) = xh(1)
         yhw(1) = yh(1)
         zhw(1) = zh(1)
         vxhw(1) = vxh(1)
         vyhw(1) = vyh(1)
         vzhw(1) = vzh(1)
         rpw(1) = rpl(1)
         rhw(1) = rhill(1)
         do i=2,njov
            j = i
            masw(j) = mjov(i)
            xhw(j) = xend(i)
            yhw(j) = yend(i)
            zhw(j) = zend(i)
            vxhw(j) = vxend(i)
            vyhw(j) = vyend(i)
            vzhw(j) = vzend(i)
            rpw(j) = rpjov(i)
            call util_hills1(mjov(1),masw(j),xhw(j),yhw(j),zhw(j),
     &            vxhw(j),vyhw(j),vzhw(j),rhw(j)) 
         enddo
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) PRIVATE(i,j)
!$OMP&   SHARED(nbod,njov,mass,xh,yh,zh,vxh,vyh,vzh,rpl,rhill,
!$OMP&          masw,xhw,yhw,zhw,vxhw,vyhw,vzhw,rpw,rhw)
         do i=2,nbod
            j = i+njov-1
            masw(j) = mass(i)
            xhw(j) = xh(i)
            yhw(j) = yh(i)
            zhw(j) = zh(i)
            vxhw(j) = vxh(i)
            vyhw(j) = vyh(i)
            vzhw(j) = vzh(i)
            rpw(j) = rpl(i)
            rhw(j) = rhill(i)
         enddo
!$OMP END PARALLEL DO

c... If it is time, output orb. elements, 
         if (t.ge.tout) then 
ccc            if (btest(iflgchk,1))  then ! bit 1 is set
               call io_write_frame_r_symba(t,ntot,nbodm,njov,masw,
     &              mjov(1),xhw,yhw,zhw,vxhw,vyhw,vzhw,outfile,
     &              iub,fopenstat)
               call io_write_mass_r(t,ntot,masw,outfile,ium,fopenstat)
ccc            endif
            if(btest(iflgchk,2))  then ! bit 2 is set
               call anal_energy_write(t,ntot,masw,j2rp2,j4rp4,
     &             xhw,yhw,zhw,vxhw,vyhw,vzhw,iue,fopenstat,eoff)
            endif
            tout = tout + dtout
         endif

c... If it is time, do a dump
         if (t.ge.tdump) then
            tfrac = (t-t0)/(tstop-t0)
            write(*,998) t,tfrac,ntot-1,nbodm-1,njov-1
 998        format(' T = ',1p1e12.5,': done = ',0pf5.3,
     &      ': Pl.tot = ',i5,': Pl.terr = ',i3,': Pl.jov = ',i3)
c            call flush()
            call io_dump_pl_symba(filedumppl,ntot,masw,xhw,yhw,zhw,
     &            vxhw,vyhw,vzhw,lclose,iflgchk,rpw,rhw,j2rp2,j4rp4)
            call io_dump_param(filedumpparam,t,tstop,dt,dtout,
     &         dtdump,iflgchk,rmin,rmax,rmaxu,qmin,lclose,outfile)
            tdump = tdump + dtdump
         endif

c... Assign planetary position and masses at the beginning of the next step
c... only jovians
         mjbeg(1) = mjov(1)
         do ipl=2,njov
            mjbeg(ipl) = mjov(ipl)
            xbeg(ipl) = xend(ipl)
            ybeg(ipl) = yend(ipl)
            zbeg(ipl) = zend(ipl)
            vxbeg(ipl) = vxend(ipl)
            vybeg(ipl) = vyend(ipl)
            vzbeg(ipl) = vzend(ipl)
         enddo

      enddo
c********** end of the big loop from time 't0' to time 'tstop'

c... Do a final dump for possible resumption later 

      call io_dump_pl_symba(filedumppl,ntot,masw,xhw,yhw,zhw,
     &          vxhw,vyhw,vzhw,lclose,iflgchk,rpw,rhw,j2rp2,j4rp4)
      call io_dump_param(filedumpparam,t,tstop,dt,dtout,
     &       dtdump,iflgchk,rmin,rmax,rmaxu,qmin,lclose,outfile)

c... get final time
c      tfin = time8()
c      write(*,*) 'Execution time',float(tfin-tini)/3600.,' h'
c      write(*,*) ctime(tfin)

      goto 666  ! diabolic trick by Morby

 2    write(*,*) 'Error or end of planet files at time=',t
      write(*,*) '  Check if time is real4 or real8 in io_read subs'
      call util_exit(1)

 666  call util_exit(0)
      close(90)
      close(91)

      end  

c*************************************************************************
c                            SYMBA7_STEP_PL_BARY.F
c*************************************************************************
c
c             Input:
c                 time          ==>  current time (real scalar)
c                 nbod          ==>  number of massive bodies (int scalar)
c                 nbodm         ==>  location of the last massie body 
c                                    (int scalar)
c                 mass          ==>  mass of bodies (real array)
c                 j2rp2,j4rp4   ==>  J2*radii_pl^2 and  J4*radii_pl^4
c                                     (real scalars)
c                 xh,yh,zh      ==>  initial position in helio coord 
c                                    (real arrays)
c                 vxh,vyh,vzh   ==>  initial velocity in helio coord 
c                                    (real arrays)
c                 vxb,vyb,vzb   ==>  initial velocity in bary coord 
c                                    (real arrays)
c                 dt            ==>  time step
c                 lclose        ==> .true. --> marge particles if they
c                                    get too close. Read in that 
c                                    distance in io_init_pl
c                                      (logical*2 scalar)
c                 rpl           ==>  physical size of a planet.
c                                    (real array)
c                 eoff          ==>  Energy offset (real scalar)
c                 rhill         ==>  size of planet's hills sphere
c                                    (real array)
c                 mtiny         ==>  Small mass  (real array)
c             Output:
c                 xh,yh,zh      ==>  final position in helio coord 
c                                       (real arrays)
c                 vxb,vyb,vzb   ==>  final velocity in bary coord 
c                                       (real arrays)
c                 rpl           ==>  Recalculated physical size of a planet.
c                                    if merger happened (real array)
c                 nbod          ==>  Recalculated number of massive bodies 
c                                    if merger happened (int scalar)
c                 mass          ==>  Recalculated mass of bodies 
c                                    if merger happened (real array)
c                 isenc         ==>  0 --> No encounter during last dt
c                                    1 --> There was encounters
c                                     (integer scalar)
c                 mergelst      ==>  list of mergers (int array)
c                 mergecnt      ==>  count of mergers (int array)
c                 iecnt         ==>  Number of encounters (int*2 array)
c                 eoff          ==>  Energy offset (real scalar)
c                 rhill         ==>  size of planet's hills sphere
c                                    (real array)
c
c Remarks: Based on symba7_step_pl - I/O bary velocitites

      subroutine symba7_step_pl_bary(time,nbod,nbodm,mass,j2rp2,
     &     j4rp4,xh,yh,zh,vxh,vyh,vzh,vxb,vyb,vzb,dt,lclose,rpl,
     &     isenc,mergelst,mergecnt,iecnt,eoff,rhill,mtiny)

      include './swift.inc'
      include './symba7.inc'

c...  Inputs Only: 
      integer nbod,nbodm
      real*8 mass(nbod),dt,time,j2rp2,j4rp4,mtiny
      logical*2 lclose 

c...  Inputs and Outputs:
      real*8 xh(nbod),yh(nbod),zh(nbod)
      real*8 vxh(nbod),vyh(nbod),vzh(nbod)
      real*8 vxb(nbod),vyb(nbod),vzb(nbod)
      real*8 rpl(nbod),eoff,rhill(NTPMAX)

c...  Outputs only
      integer isenc
      integer iecnt(NTPMAX),ielev(NTPMAX)
      integer mergelst(2,NTPMAX),mergecnt

c...  Internals
      integer i,j,ieflg,irec
      logical*1 svdotr            ! Not used in the routine
      integer ielst(2,NENMAX),ielc

c----
c...  Executable code 

      do i=2,nbod
         iecnt(i) = 0
         ielev(i) = -1
      enddo
      isenc = 0

c...  check for encounters
      irec = 0
      ielc = 0
      do j=2,nbodm
!$OMP       PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) 
!$OMP&      PRIVATE(i,ieflg,svdotr)
!$OMP&      SHARED(j,nbod,rhill,mass,xh,yh,zh,vxh,vyh,vzh,dt,
!$OMP&             irec,isenc,iecnt,ielev,ielc,ielst)
         do i=j+1,nbod
            ieflg = 0
            call symba7_chk(rhill,nbod,i,j,mass,xh,yh,zh,vxh,
     &           vyh,vzh,dt,irec,ieflg,svdotr)
            if(ieflg.ne.0) then
!$OMP          CRITICAL
               isenc = 1
               iecnt(i) = iecnt(i) + 1
               iecnt(j) = iecnt(j) + 1
               ielev(i) = 0
               ielev(j) = 0
               ielc = ielc + 1
               if(ielc.gt.NENMAX) then
                  write(*,*) 'ERROR: encounter matrix is filled.'
                  write(*,*) 'STOPPING'
                  call util_exit(1)
               endif
               ielst(1,ielc) = i
               ielst(2,ielc) = j
!$OMP          END CRITICAL
            endif
         enddo
!$OMP    END PARALLEL DO
      enddo

c...  do a step
      if(isenc.eq.0) then
         call symba7_step_helio(nbod,nbodm,mass,j2rp2,
     &     j4rp4,xh,yh,zh,vxb,vyb,vzb,dt)
         mergecnt=0
      else
         call symba7_step_interp(time,iecnt,ielev,nbod,
     &     nbodm,mass,rhill,j2rp2,j4rp4,lclose,rpl,xh,yh,zh,
     &     vxb,vyb,vzb,dt,mergelst,mergecnt,eoff,ielc,ielst,
     &     mtiny)  
      endif

      return
      end 

c*************************************************************************
c                            SYMBA7_STEP_HELIO.F
c*************************************************************************
c This subroutine takes a step in helio coord.  
c Does a KICK than a DRIFT than a KICK.
c ONLY DOES MASSIVE PARTICLES
c
c             Input:
c                 nbod          ==>  number of massive bodies (int scalar)
c                 mass          ==>  mass of bodies (real array)
c                 j2rp2,j4rp4   ==>  J2*radii_pl^2 and  J4*radii_pl^4
c                                     (real scalars)
c                 xh,yh,zh      ==>  initial position in helio coord 
c                                    (real arrays)
c                 vxb,vyb,vzb   ==>  initial velocity in bary coord 
c                                    (real arrays)
c                 dt            ==>  time step
c             Output:
c                 xh,yh,zh      ==>  final position in helio coord 
c                                       (real arrays)
c                 vxb,vyb,vzb   ==>  final velocity in bary coord 
c                                       (real arrays)
c
c Remarks: Based on helio_step_pl.f but does not pass the intermediate
c          positions and velocities back for the TP to use.
c          Input and output velocities are barycentric

      subroutine symba7_step_helio(nbod,nbodm,mass,j2rp2,
     &     j4rp4,xh,yh,zh,vxb,vyb,vzb,dt)

      include './swift.inc'

c...  Inputs Only: 
      integer nbod,nbodm
      real*8 mass(nbod),dt,j2rp2,j4rp4

c...  Inputs and Outputs:
      real*8 xh(nbod),yh(nbod),zh(nbod)
      real*8 vxb(nbod),vyb(nbod),vzb(nbod)

c...  Internals:
      integer i1stloc
      real*8 dth
      real*8 axh(NTPMAX),ayh(NTPMAX),azh(NTPMAX)
      real*8 ptxb,ptyb,ptzb            ! Not used here
      real*8 ptxe,ptye,ptze

c----
c...  Executable code 

      dth = 0.5d0*dt
      i1stloc = 0    ! always compute the accels at the begining and end of step

c...  Do the linear drift due to momentum of the Sun
      call helio_lindrift(nbod,mass,vxb,vyb,vzb,dth,
     &     xh,yh,zh,ptxb,ptyb,ptzb)

c...  Get the accelerations in helio frame. if frist time step
      call symba7_helio_getacch(i1stloc,nbod,nbodm,mass,j2rp2,j4rp4,
     &     xh,yh,zh,axh,ayh,azh)

c...  Apply a heliocentric kick for a half dt 
      call kickvh(nbod,vxb,vyb,vzb,axh,ayh,azh,dth)

c..   Drift in helio coords for the full step 
      call helio_drift(nbod,mass,xh,yh,zh,vxb,vyb,vzb,dt)

c...  Get the accelerations in helio frame.
      call symba7_helio_getacch(i1stloc,nbod,nbodm,mass,j2rp2,j4rp4,
     &     xh,yh,zh,axh,ayh,azh)

c...  Apply a heliocentric kick for a half dt 
      call kickvh(nbod,vxb,vyb,vzb,axh,ayh,azh,dth)

c...  Do the linear drift due to momentum of the Sun
      call helio_lindrift(nbod,mass,vxb,vyb,vzb,dth,
     &     xh,yh,zh,ptxe,ptye,ptze)

      return
      end   

c*************************************************************************
c                            SYMBA7_STEP_INTERP.F
c*************************************************************************
c
c             Input:
c                 time          ==> Current time (real scalar)
c                 iecnt         ==>  The number of objects that each planet 
c                                    is encountering (int*2 array)
c                 ielev         ==>  The level that this particle should go
c                                             (int*2 array)
c                 nbod          ==>  number of massive bodies (int scalar)
c                 nbodm         ==>  Location of last massive body(int scalar)
c                 mass          ==>  mass of bodies (real array)
c                 rhill         ==>  Radius of hill sphere (real array)
c                 j2rp2,j4rp4   ==>  J2*radii_pl^2 and  J4*radii_pl^4
c                                     (real scalars)
c                 xh,yh,zh      ==>  initial position in helio coord 
c                                    (real arrays)
c                 vxb,vyb,vzb   ==>  initial velocity in bary coord 
c                                    (real arrays)
c                 dt            ==>  time step
c                 lclose        ==> .true. --> marge particles if they
c                                    get too close. Read in that 
c                                    distance in io_init_pl
c                                      (logical*2 scalar)
c                 rpl           ==>  physical size of a planet.
c                                    (real array)
c                 eoff          ==>  Energy offset (real scalar)
c                ielc           ==>  number of encounters (integer*2 scalar)
c                ielst          ==>  list of ecnounters (2D integer*2 array)
c                mtiny          ==>  Small mass  (real array)
c             Output:
c                 xh,yh,zh      ==>  final position in helio coord 
c                                       (real arrays)
c                 vxb,vyb,vzb   ==>  final velocity in bary coord 
c                                       (real arrays)
c                 rpl           ==>  Recalculated physical size of a planet.
c                                    if merger happened (real array)
c                 nbod          ==>  Recalculated number of massive bodies 
c                                    if merger happened (int scalar)
c                 nbodm         ==>  Location of last massive body(int scalar)
c                 mass          ==>  Recalculated mass of bodies 
c                                    if merger happened (real array)
c                 mergelst      ==>  list of mergers (int array)
c                 mergecnt      ==>  count of mergers (int array)
c                 eoff          ==>  Energy offset (real scalar)
c
c Remarks: Input and output velocities are barycentric

      subroutine symba7_step_interp(time,iecnt,ielev,nbod,
     &     nbodm,mass,rhill,j2rp2,j4rp4,lclose,rpl,xh,yh,zh,vxb,
     &     vyb,vzb,dt,mergelst,mergecnt,eoff,ielc,ielst,mtiny)

      include './swift.inc'
      include './symba7.inc'

c...  Inputs Only: 
      real*8 mass(NTPMAX),dt,j2rp2,j4rp4,time,mtiny
      integer iecnt(NTPMAX),ielev(NTPMAX)
      logical*2 lclose 
      integer ielst(2,NENMAX),ielc

c...  Inputs and Outputs:
      integer nbod,nbodm
      real*8 xh(nbod),yh(nbod),zh(nbod)
      real*8 vxb(nbod),vyb(nbod),vzb(nbod)
      real*8 rpl(nbod),eoff
      real*8 rhill(nbod)

c...  Outputs
      integer mergelst(2,NTPMAX),mergecnt

c...  Internals:
      integer irec,ilevl(NTPMAX),i 
      real*8 dth
      real*8 axh(NTPMAX),ayh(NTPMAX),azh(NTPMAX)
      real*8 ptxb,ptyb,ptzb            ! Not used here
      real*8 ptxe,ptye,ptze
      logical*1 svdotr(NENMAX)  ! Used by symba_step_recur

      save axh,ayh,azh     ! Note this !!

c----
c...  Executable code 

      dth = 0.5d0*dt

c...  Do the linear drift due to momentum of the Sun
      call helio_lindrift(nbod,mass,vxb,vyb,vzb,dth,
     &     xh,yh,zh,ptxb,ptyb,ptzb)

c...  Get the accelerations in helio frame. For each object
c...     only include those guys that it is not encountering with. 
      call symba7_getacch(nbod,nbodm,mass,j2rp2,j4rp4,
     &     xh,yh,zh,axh,ayh,azh,mtiny,ielc,ielst)

c...  Apply a heliocentric kick for a half dt 
      call kickvh(nbod,vxb,vyb,vzb,axh,ayh,azh,dth)

c..   Do a recursion step for full dt
      irec = -1
      call symba7_helio_drift(nbod,ielev,irec,mass,xh,
     &     yh,zh,vxb,vyb,vzb,dt)
      irec = 0
      do i=2,nbod
         ilevl(i) = 0
      enddo
      mergecnt = 0
      call symba7_step_recur(time,nbod,nbodm,mass,irec,ilevl,
     &     iecnt,ielev,rhill,xh,yh,zh,vxb,vyb,vzb,lclose,
     &     rpl,mergelst,mergecnt,dt,eoff,svdotr,ielc,ielst)

c...  Get the accelerations in helio frame. For each object
c...     only include those guys that it is not encountering with. 
      call symba7_getacch(nbod,nbodm,mass,j2rp2,j4rp4,
     &     xh,yh,zh,axh,ayh,azh,mtiny,ielc,ielst)

c...  Apply a heliocentric kick for a half dt 
      call kickvh(nbod,vxb,vyb,vzb,axh,ayh,azh,dth)

c...  Do the linear drift due to momentum of the Sun
      call helio_lindrift(nbod,mass,vxb,vyb,vzb,dth,
     &     xh,yh,zh,ptxe,ptye,ptze)

      return
      end   

c***********************************************************************
c                          COORD_VH2B_JOV.F
c***********************************************************************
c     PURPOSE: Converts from Heliocentric to Barycentric coords. 
c              Velocity only - Jovians only
c     ARGUMENTS:  Input is 
c                    ibar = 0 => input for terrestrial is helio
c                         = 1 => input for terrestrial is bary
c                    nbod => number of terrestrials (integer)
c                    mass(*) =>  masses of terrestrials (real array)
c                    njov => number of jovians (integer)
c                    massj(*) =>  masses of jovians (real array)
c                    vx(*),vy(*),vz(*) => Terrestrials velocities (real array)
c                                         helio or bary depending on ibar
c                    vxhj(*),vyhj(*),vzhj(*) => Jovians helio velocities
c                                             (real array)
c                 Returned are
c                    vxbj(*),vybj(*),vzbj(*) => Jovian bary velocities
c                                            (real array)
c Remarks: Adapted from coord_vh2b
c          PARALELIZED

      subroutine coord_vh2b_jov(ibar,nbod,mass,vx,vy,vz,
     &      njov,massj,vxhj,vyhj,vzhj,vxbj,vybj,vzbj)

      include './swift.inc'

c...  Inputs: 
      integer nbod,njov,ibar
      real*8 mass(nbod),massj(njov)
      real*8 vx(nbod),vy(nbod),vz(nbod)
      real*8 vxhj(njov),vyhj(njov),vzhj(njov)

c...  Outputs:
      real*8 vxbj(njov),vybj(njov),vzbj(njov)

c...  Internals:
      integer n
      real*8 msys,msysj
      real*8 vxtmp,vytmp,vztmp

c----
c...  Executable code 

c...  First add up the jovians
      msysj = mass(1) + massj(2)
      vxbj(1) = massj(2)*vxhj(2)
      vybj(1) = massj(2)*vyhj(2)
      vzbj(1) = massj(2)*vzhj(2)
      do n=3,njov
         msysj = msysj + massj(n)
         vxbj(1) = vxbj(1) + massj(n)*vxhj(n)
         vybj(1) = vybj(1) + massj(n)*vyhj(n)
         vzbj(1) = vzbj(1) + massj(n)*vzhj(n)
      enddo

c...  Then add up the terrestrials
      msys = 0.0d0
      vxtmp = 0.0d0
      vytmp = 0.0d0
      vztmp = 0.0d0
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) 
!$OMP&    PRIVATE(n) SHARED(nbod,mass,vx,vy,vz) 
!$OMP&    REDUCTION(+:msys,vxtmp,vytmp,vztmp) 
      do n=2,nbod
         msys = msys + mass(n)
         vxtmp = vxtmp + mass(n)*vx(n)
         vytmp = vytmp + mass(n)*vy(n)
         vztmp = vztmp + mass(n)*vz(n)
      enddo
!$OMP END PARALLEL DO
      msys = msys + msysj
      vxbj(1) = vxbj(1) + vxtmp 
      vybj(1) = vybj(1) + vytmp 
      vzbj(1) = vzbj(1) + vztmp 


c...  Compute the barycentric velocity of the Sun
      if (ibar.eq.0) then
         vxbj(1) = -vxbj(1)/msys
         vybj(1) = -vybj(1)/msys
         vzbj(1) = -vzbj(1)/msys
      else
         vxbj(1) = -vxbj(1)/msysj
         vybj(1) = -vybj(1)/msysj
         vzbj(1) = -vzbj(1)/msysj
      endif

c...  Barycentric velocities of jovians
      do n=2,njov
        vxbj(n) = vxhj(n) + vxbj(1)
        vybj(n) = vyhj(n) + vybj(1)
        vzbj(n) = vzhj(n) + vzbj(1)
      enddo

      return
      end     ! coord_vh2b_jov

c***********************************************************************
c                          COORD_VH2B_TER.F
c***********************************************************************
c     PURPOSE: Converts from Helio to Bary coords.
c              Velocity only - Terrestrials only
c     ARGUMENTS:  Input is 
c                    nbod => number of terrestrials (integer)
c                    njov => number of jovians (integer)
c                    vxh(*),vyh(*),vzh(*) => Terrestrials helio velocities
c                                             (real array)
c                    vxbj(*),vybj(*),vzbj(*) => Jovians bary velocities
c                                             (real array)
c                 Returned are
c                    vxb(*),vyb(*),vzb(*) => Terrestrials bary velocities
c                                            (real array)
c Remarks: Adapted from coord_vh2b
c          PARALELIZED

      subroutine coord_vh2b_ter(nbod,vxh,vyh,vzh,njov,
     &             vxbj,vybj,vzbj,vxb,vyb,vzb)

      include './swift.inc'

c...  Inputs: 
      integer nbod,njov
      real*8 vxh(nbod),vyh(nbod),vzh(nbod)
      real*8 vxbj(njov),vybj(njov),vzbj(njov)

c...  Outputs:
      real*8 vxb(nbod),vyb(nbod),vzb(nbod)

c...  Internals:
      integer n

c----
c...  Executable code 

!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) PRIVATE(n)
!$OMP&   SHARED(nbod,vxb,vyb,vzb,vxh,vyh,vzh,vxbj,vybj,vzbj)
      do n=2,nbod
         vxb(n) = vxh(n) + vxbj(1)
         vyb(n) = vyh(n) + vybj(1)
         vzb(n) = vzh(n) + vzbj(1)
      enddo
!$OMP END PARALLEL DO

      return
      end     ! coord_vh2b_ter

c***********************************************************************
c                          COORD_VB2H_TER.F
c***********************************************************************
c     PURPOSE: Converts from Barycentric to Helio coords.
c              Velocity only - Terrestrials only
c     ARGUMENTS:  Input is 
c                    nbod => number of terrestrials (integer)
c                    mass(*) =>  masses of terrestrials (real array)
c                    njov => number of jovians (integer)
c                    massj(*) =>  masses of jovians (real array)
c                    vxb(*),vyb(*),vzb(*) => Terrestrials bary velocities (real array)
c                    vxbj(*),vybj(*),vzbj(*) => Jovians bary velocities
c                                             (real array)
c                 Returned are
c                    vxh(*),vyh(*),vzh(*) => Terrestrials helio velocities
c                                            (real array)

      subroutine coord_vb2h_ter(nbod,mass,vxb,vyb,vzb,njov,massj,
     &             vxbj,vybj,vzbj,vxh,vyh,vzh)

      include './swift.inc'

c...  Inputs: 
      integer nbod,njov
      real*8 mass(nbod),massj(njov)
      real*8 vxb(nbod),vyb(nbod),vzb(nbod)
      real*8 vxbj(njov),vybj(njov),vzbj(njov)

c...  Outputs:
      real*8 vxh(nbod),vyh(nbod),vzh(nbod)

c...  Internals:
      integer i
      real*8 vxtmp,vytmp,vztmp

c----
c...  Executable code 

c...  First add up the jovians
      vxbj(1) = -massj(2)*vxbj(2)
      vybj(1) = -massj(2)*vybj(2)
      vzbj(1) = -massj(2)*vzbj(2)
      do i=3,njov
         vxbj(1) = vxbj(1) - massj(i)*vxbj(i)
         vybj(1) = vybj(1) - massj(i)*vybj(i)
         vzbj(1) = vzbj(1) - massj(i)*vzbj(i)
      enddo

c...  Then add up the terrestrials
      vxtmp = 0.0d0
      vytmp = 0.0d0
      vztmp = 0.0d0
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) PRIVATE(i)
!$OMP&   SHARED(nbod,mass,vxb,vyb,vzb)
!$OMP&   REDUCTION(+:vxtmp,vytmp,vztmp)
      do i=2,nbod
         vxtmp = vxtmp - mass(i)*vxb(i)
         vytmp = vytmp - mass(i)*vyb(i)
         vztmp = vztmp - mass(i)*vzb(i)
      enddo
!$OMP END PARALLEL DO
      vxbj(1) = vxbj(1) + vxtmp 
      vybj(1) = vybj(1) + vytmp 
      vzbj(1) = vzbj(1) + vztmp 

c...  Helio velocities of terrestrials
      vxbj(1) = vxbj(1)/mass(1)
      vybj(1) = vybj(1)/mass(1)
      vzbj(1) = vzbj(1)/mass(1)

!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) PRIVATE(i)
!$OMP&   SHARED(nbod,vxb,vyb,vzb,vxh,vyh,vzh,vxbj,vybj,vzbj)
      do i=2,nbod
         vxh(i) = vxb(i) - vxbj(1)  
         vyh(i) = vyb(i) - vybj(1)  
         vzh(i) = vzb(i) - vzbj(1)  
      enddo
!$OMP END PARALLEL DO

      return
      end     ! coord_vb2h_ter

c*************************************************************************
c                            HELIO_LINDRIFT_JOV.F
c*************************************************************************
c This subroutine takes a linear drift due to mometum of Sun
c
c             Input:
c                 nbod          ==>  number of massive bodies (int scalar)
c                 mass          ==>  mass of bodies (real array)
c                 vxb,vyb,vzb   ==>  velocity in bary coord 
c                                    (real arrays)
c                 dt            ==>  time step
c                 xh,yh,zh      ==>  initial position in helio coord 
c                                       (real arrays)
c             Output:
c                 xh,yh,zh      ==>  final position in helio coord 
c                                       (real arrays)
c
c Remarks: Bases on Martin's code h2.f
c Authors:  Hal Levison 
c Date:    11/14/96
c Last revision: 1/8/97

      subroutine helio_lindrift_jov(nbod,mass,njov,massj,
     &     vxbj,vybj,vzbj,dt,xh,yh,zh)

      include './swift.inc'

c...  Inputs Only: 
      integer nbod,njov
      real*8 mass(nbod),massj(njov),dt
      real*8 vxbj(njov),vybj(njov),vzbj(njov)

c...  Inputs and Outputs:
      real*8 xh(nbod),yh(nbod),zh(nbod)

c...  Internals:
      integer n
      real*8 ptx,pty,ptz

c----
c...  Executable code 

c...  First add the terrestrials
      ptx = 0.0d0
      pty = 0.0d0
      ptz = 0.0d0

c...  Then add the jovians
      do n=2,njov
         ptx = ptx + massj(n)*vxbj(n)
         pty = pty + massj(n)*vybj(n)
         ptz = ptz + massj(n)*vzbj(n)
      enddo

      ptx = ptx/mass(1)
      pty = pty/mass(1)
      ptz = ptz/mass(1)

c...  Apply a drift on terrestrials only
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) PRIVATE(n)
!$OMP&   SHARED(nbod,xh,yh,zh,ptx,pty,ptz,dt)
      do n=2,nbod
         xh(n) = xh(n) + ptx*dt
         yh(n) = yh(n) + pty*dt
         zh(n) = zh(n) + ptz*dt
      enddo
!$OMP END PARALLEL DO

      return
      end       ! helio_lindrift_jov

c*************************************************************************
c                        HELIO_GETACCH_JOV.F
c*************************************************************************
c This subroutine calculates the 3rd term acceleration on the test particles
c in the HELIOCENTRIC frame. This term is the direct cross terms
c             Input:
c                 nbod          ==>  number of massive bodies (int scalor)
c                 ntp           ==>  number of test particles (int scalor)
c                 mass          ==>  mass of massive bodies (real array)
c                 xh,yh,zh      ==>  massive position in heliocentric coord 
c                                   (real arrays)
c                 xht,yht,zht   ==>  tp position in heliocentric coord 
c                                   (real arrays)
c                  istat       ==>  status of the test paricles
c                                      (integer array)
c                                      istat(i) = 0 ==> active:  = 1 not
c                                    NOTE: it is really a 2d array but 
c                                          we only use the 1st row
c             Output:
c                 axht,ayht,azht ==>  3rd term of acceleration in helio coord 
c                                     (real arrays)
c
c Author:  Hal Levison  
c Date:    2/18/93
c Last revision: 

      subroutine helio_getacch_jov(nbod,mass,xh,yh,zh,ntp,xht,yht,zht,
     &                     axht,ayht,azht)

      include './swift.inc'

c...  Inputs: 
      integer nbod,ntp
      real*8 mass(nbod),xh(nbod),yh(nbod),zh(nbod)
      real*8 xht(ntp),yht(ntp),zht(ntp)

c...  Outputs:
      real*8 axht(NTPMAX),ayht(NTPMAX),azht(NTPMAX)

c...  Internals:
      integer i,j
      real*8 dx,dy,dz,rji2,fac,irij3

c------
c...  Executable code

!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) 
!$OMP&   PRIVATE(i) SHARED(ntp,axht,ayht,azht)
      do i=1,ntp
         axht(i) = 0.0d0
         ayht(i) = 0.0d0
         azht(i) = 0.0d0
      enddo
!$OMP END PARALLEL DO

      do i=2,nbod
!$OMP    PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) 
!$OMP&      PRIVATE(j,dx,dy,dz,rji2,irij3,fac) 
!$OMP&      SHARED(i,ntp,mass,xht,yht,zht,xh,yh,zh,axht,ayht,azht)
         do j=2,ntp
            dx = xht(j) - xh(i)
            dy = yht(j) - yh(i)
            dz = zht(j) - zh(i)
            rji2 = dx*dx + dy*dy + dz*dz

            irij3 = 1.0d0/(rji2*sqrt(rji2))
            fac = mass(i)*irij3

            axht(j) = axht(j) - fac*dx
            ayht(j) = ayht(j) - fac*dy
            azht(j) = azht(j) - fac*dz
         enddo
!$OMP    END PARALLEL DO
      enddo

      return
      end     ! helio_getacch_jov

c*************************************************************************
c                            RMVS3_INTERP_MOD.F
c*************************************************************************
c This subroutine interpolates between two kepler orbits.
c For outer region only
c
c             Input:
c                 nbod                ==>  number of massive bodies 
c                                          (int scalar)
c                 xbeg,ybeg,zbeg      ==>  initial planet position in helio 
c                                            (real arrays)
c                 vxbeg,vybeg,vzbeg   ==>  initial planet vel in helio 
c                                            (real arrays)
c                 xend,yend,zend      ==>  final planet position in helio 
c                                            (real arrays)
c                 vxend,vyend,vzend   ==>  final planet position in helio 
c                                            (real arrays)
c                 dt                   ==>  time step (reswift_rmvs3_hybrid_tr8.fal sclar)
c                 msun                 ==>  mass of sun (real sclar)
c                 nt                   ==>  the number of intermediate steps
c                                           (integer scalar)
c             Output:
c                 xtmp,ytmp,ztmp      ==>  position of planet wrt time 
c                                          for outer region
c                                            (2d real arrays)
c                 vxtmp,vytmp,vztmp   ==>  velocoty of planet wrt time 
c                                          for outer region
c                                            (2d real arrays)
c
c Remarks: Based on rmvs2_interp_o 

      subroutine rmvs3_interp_mod(nbod,xbeg,ybeg,zbeg,vxbeg,vybeg,
     &     vzbeg,xend,yend,zend,vxend,vyend,vzend,dt,massb,masse,nt,
     &     xtmp,ytmp,ztmp,vxtmp,vytmp,vztmp)

      include './swift.inc'
      include './hybrid.inc'

c...  Inputs Only: 
      integer nbod,nt
      real*8 dt
      real*8 massb(nbod),masse(nbod)
      real*8 xbeg(nbod),ybeg(nbod),zbeg(nbod)
      real*8 vxbeg(nbod),vybeg(nbod),vzbeg(nbod)
      real*8 xend(nbod),yend(nbod),zend(nbod)
      real*8 vxend(nbod),vyend(nbod),vzend(nbod)

c...  Outputs:
      real*8 xtmp(NPLMAX,0:NINTERP),ytmp(NPLMAX,0:NINTERP)
      real*8 vxtmp(NPLMAX,0:NINTERP),vytmp(NPLMAX,0:NINTERP)
      real*8 ztmp(NPLMAX,0:NINTERP),vztmp(NPLMAX,0:NINTERP)

c...  Internals
      integer i,iflg,ib
      real*8 xc2(NPLMAX),yc2(NPLMAX),zc2(NPLMAX)
      real*8 vxc2(NPLMAX),vyc2(NPLMAX),vzc2(NPLMAX)
      real*8 xc1(NPLMAX),yc1(NPLMAX),zc1(NPLMAX)
      real*8 vxc1(NPLMAX),vyc1(NPLMAX),vzc1(NPLMAX)
      real*8 dti,dtb,frac,onemf,gm

c...  Executable code 
      dti = dt/float(nt)
      dtb = -1.0d0*dt

c...  move the end positions to beginning
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) 
!$OMP&   PRIVATE(i,gm,iflg) 
!$OMP&   SHARED(nbod,masse,dtb,
!$OMP&          xend,yend,zend,vxend,vyend,vzend,
!$OMP&          xc2,yc2,zc2,vxc2,vyc2,vzc2)
      do i=2,nbod
         xc2(i) = xend(i)
         yc2(i) = yend(i)
         zc2(i) = zend(i)
         vxc2(i) = vxend(i)
         vyc2(i) = vyend(i)
         vzc2(i) = vzend(i)

c...  uses the mass at the end of the step
         gm = masse(1)+masse(i)
         call drift_one(gm,xc2(i),yc2(i),zc2(i),
     &             vxc2(i),vyc2(i),vzc2(i),dtb,iflg)
         if (iflg.ne.0) then
!$OMP       CRITICAL
            write(*,*) ' Planet',i,' is lost in rmvs3_interp_mod !'
            write(*,*) masse(1),dtb
            write(*,*) xc2(i),yc2(i),zc2(i)
            write(*,*) vxc2(i),vyc2(i),vzc2(i)
            write(*,*) ' STOPPING '
            call util_exit(1)
!$OMP       END CRITICAL
         endif
      enddo
!$OMP END PARALLEL DO

c...  compute the intermediate steps
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) 
!$OMP&   PRIVATE(i,ib,gm,iflg,frac,onemf) 
!$OMP&   SHARED(nbod,xtmp,ytmp,ztmp,vxtmp,vytmp,vztmp,
!$OMP&          xbeg,ybeg,zbeg,vxbeg,vybeg,vzbeg,nt,
!$OMP&          massb,xc1,yc1,zc1,vxc1,vyc1,vzc1,dti,
!$OMP&          masse,xc2,yc2,zc2,vxc2,vyc2,vzc2,dtb)
      do i=2,nbod
         xc1(i) = xbeg(i)
         yc1(i) = ybeg(i)
         zc1(i) = zbeg(i)
         vxc1(i) = vxbeg(i)
         vyc1(i) = vybeg(i)
         vzc1(i) = vzbeg(i)

         xtmp(i,0) = xbeg(i)
         ytmp(i,0) = ybeg(i)
         ztmp(i,0) = zbeg(i)
         vxtmp(i,0) = vxbeg(i)
         vytmp(i,0) = vybeg(i)
         vztmp(i,0) = vzbeg(i)

         do ib = 1,nt

c...  use the mass at the begin of the step
            gm = massb(1)+massb(i)
            call drift_one(gm,xc1(i),yc1(i),zc1(i),
     &           vxc1(i),vyc1(i),vzc1(i),dti,iflg)
            if (iflg.ne.0) then
!$OMP          CRITICAL
               write(*,*) ' Planet',i,' is lost in rmvs3_interp_mod !'
               write(*,*) massb(1),dtb
               write(*,*) xc1(i),yc1(i),zc1(i)
               write(*,*) vxc1(i),vyc1(i),vzc1(i)
               write(*,*) ' STOPPING '
               call util_exit(1)
!$OMP          END CRITICAL
            endif

c...  use the mass at the end of the step
            gm = masse(1)+masse(i)
            call drift_one(gm,xc2(i),yc2(i),zc2(i),
     &           vxc2(i),vyc2(i),vzc2(i),dti,iflg)
            if (iflg.ne.0) then
!$OMP          CRITICAL
               write(*,*) ' Planet',i,' is lost in rmvs3_interp_mod !'
               write(*,*) masse(1),dtb
               write(*,*) xc2(i),yc2(i),zc2(i)
               write(*,*) vxc2(i),vyc2(i),vzc2(i)
               write(*,*) ' STOPPING '
               call util_exit(1)
!$OMP          END CRITICAL
            endif

c...  perform the interpolation
            frac = float(ib)/float(nt)
            onemf = 1.0d0 - frac
            xtmp(i,ib) = onemf*xc1(i) + frac*xc2(i)
            ytmp(i,ib) = onemf*yc1(i) + frac*yc2(i)
            ztmp(i,ib) = onemf*zc1(i) + frac*zc2(i)
            vxtmp(i,ib) = onemf*vxc1(i) + frac*vxc2(i)
            vytmp(i,ib) = onemf*vyc1(i) + frac*vyc2(i)
            vztmp(i,ib) = onemf*vzc1(i) + frac*vzc2(i)
         enddo
      enddo
!$OMP END PARALLEL DO

c...  put zeros in position 1 (Sun)
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) 
!$OMP&   PRIVATE(ib) 
!$OMP&   SHARED(nt,xtmp,ytmp,ztmp,vxtmp,vytmp,vztmp)
      do ib = 0,nt
         xtmp(1,ib) = 0.0d0
         ytmp(1,ib) = 0.0d0
         ztmp(1,ib) = 0.0d0
         vxtmp(1,ib) = 0.0d0
         vytmp(1,ib) = 0.0d0
         vztmp(1,ib) = 0.0d0
      enddo
!$OMP END PARALLEL DO

      return
      end      ! rmvs3_interp_mod.f

c*************************************************************************
c                            IO_WRITE_FRAME_R_SYMBA
c*************************************************************************
c write out a whole frame to an real*4 binary file.
c both massive and test particles
c
c             Input:
c                 time          ==>  current time (real scalar)
c                 nbod          ==>  number of massive bodies (int scalar)
c                 ntp            ==>  number of massive bodies (int scalar)
c                 mass          ==>  mass of bodies (real array)
c                 xh,yh,zh      ==>  current position in helio coord
c                                    (real arrays)
c                 vxh,vyh,vzh   ==>  current velocity in helio coord
c                                    (real arrays)
c                 xht,yht,zht    ==>  current part position in helio coord
c                                      (real arrays)
c                 vxht,vyht,vzht ==>  current velocity in helio coord
c                                        (real arrays)
c                 istat           ==>  status of the test paricles
c                                      (2d integer array)
c                                      istat(i,1) = 0 ==> active:  = 1 not
c                                      istat(i,2) = -1 ==> Danby did not work
c                 oname           ==> output file name (character string)
c                 iu              ==> unit number to write to
c                 fopenstat       ==>  The status flag for the open
c                                      statements of the output files.
c                                          (character*80)
c
c
c Remarks: Based on io_write_frame
c Authors:  Hal Levison
c Date:    2/22/94
c Last revision:

      subroutine io_write_frame_r_symba(time,nbod,nter,njov,mass,
     &           msun,xh,yh,zh,vxh,vyh,vzh,oname,iu,fopenstat)

      include './swift.inc'
      include './io.inc'

c...  Inputs:
      integer nbod,nter,njov,iu
      real*8 mass(nbod),time,msun
      real*8 xh(nbod),yh(nbod),zh(nbod)
      real*8 vxh(nbod),vyh(nbod),vzh(nbod)
      character*80 oname,fopenstat

c...  Internals
      integer i,id
      integer ialpha,ierr
      real*8 a,e,inc,capom,omega,capm
      real*8 gm
      integer i1st    ! =0 first time through; =1 after
      data i1st/0/
      save i1st
      
c      real*4 time4
c      integer*2 nbod2,nter2,njov2

c----
c...  Executable code

c...  if first time through open file
      if(i1st.eq.0) then
         call io_open(iu,oname,fopenstat,'UNFORMATTED',ierr)
         if(ierr.ne.0) then
           write(*,*) ' SWIFT ERROR: in io_write_frame: '
           write(*,*) '     Could not open binary output file:'
           call util_exit(1)
         endif
         i1st = 1
      else
        call io_open(iu,oname,'append','UNFORMATTED',ierr)
      endif

c... write header (note the three nbod numbers)
c      nbod2 = nbod
c      nter2 = nter
c      njov2 = njov
c      time4 = time
c      write(iu) time4,nbod2,nter2,njov2
      write(iu)time,nbod,njov	

c...  write out planets
c...  first terrestrials massive
      do i=2,nter
         gm = mass(1)+mass(i)
         id = i
         call orbel_xv2el(xh(i),yh(i),zh(i),vxh(i),vyh(i),vzh(i),gm,
     &          ialpha,a,e,inc,capom,omega,capm)
         call io_write_line_r(iu,id,a,e,inc,capom,omega,capm)
      enddo
  
c...  then jovians (note different msun and index)
      do i=nter+1,nter+njov-1
         gm = msun+mass(i)
         id = i
         call orbel_xv2el(xh(i),yh(i),zh(i),vxh(i),vyh(i),vzh(i),gm,
     &          ialpha,a,e,inc,capom,omega,capm)
         call io_write_line_r(iu,id,a,e,inc,capom,omega,capm)
      enddo

c... then terrestrials planetesimals
      do i=nter+njov,nbod
         gm = mass(1)+mass(i)
         id = i
         call orbel_xv2el(xh(i),yh(i),zh(i),vxh(i),vyh(i),vzh(i),gm,
     &          ialpha,a,e,inc,capom,omega,capm)
         call io_write_line_r(iu,id,a,e,inc,capom,omega,capm)
      enddo

      close(iu)
      return
      end      

c*************************************************************************
c                        GETACCH_GR.F
c*************************************************************************
c This subroutine calculates the acceleration on the planets due to
c relativistic effects
c             Input:
c                  nbod        ==>  number of massive bodies (int scalor)
c                  mass        ==>  mass of bodies (real array)
c                  xh,yh,zh    ==>  massive part position in helio coord 
c                                     (real arrays)
c                  vxh,vyh,vzh ==>  massive part velocity in helio coord 
c                                     (real arrays)
c             Output:
c                  axh,ayh,azh ==>  massive part acceleration in helio coord 
c                                   (real arrays)
c
c Author:  F. Roig
c Date:    20/11/15

      subroutine getacch_gr(nbod,mass,xh,yh,zh,vxh,vyh,vzh,
     &     axh,ayh,azh,nmax,gmc2)

      include './swift.inc'

c...  Inputs: 
      integer nbod,nmax
      real*8 mass(nbod),xh(nbod),yh(nbod),zh(nbod)
      real*8 vxh(nbod),vyh(nbod),vzh(nbod)
      real*8 gmc2

c...  Outputs:
      real*8 axh(NPLMAX),ayh(NPLMAX),azh(NPLMAX)
                
c...  Internals:
      integer i
      real*8 ir3h(NPLMAX),irh(NPLMAX)
      real*8 rv,v2,fac1,fac2
c      integer ialpha
c      real*8 gm,a,e,q
      
c----
c...  Executable code 

c...  get thr r^-3's  for the planets
      call getacch_ir3(nmax,2,xh,yh,zh,ir3h,irh)

c...  add them all together

!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) 
!$OMP&   PRIVATE(i,rv,v2,fac1,fac2) 
!$OMP&   SHARED(nbod,mass,xh,yh,zh,vxh,vyh,vzh,ir3h,irh,
!$OMP&          axh,ayh,azh,nmax,gmc2)
      do i=2,nbod
         axh(i) = 0.0d0
         ayh(i) = 0.0d0
         azh(i) = 0.0d0
         
         if (i.le.nmax) then

c... formula of Quinn et al 1991
         rv = xh(i)*vxh(i) + yh(i)*vyh(i) + zh(i)*vzh(i)
         v2 = vxh(i)*vxh(i) + vyh(i)*vyh(i) + vzh(i)*vzh(i)
         fac1 = gmc2*ir3h(i)
         fac2 = 4.d0*mass(1)*irh(i) - v2
         axh(i) = fac1*(fac2*xh(i) + 4.d0*rv*vxh(i))
         ayh(i) = fac1*(fac2*yh(i) + 4.d0*rv*vyh(i))
         azh(i) = fac1*(fac2*zh(i) + 4.d0*rv*vzh(i))

c... formula of Nobili & Roxburgh 1986
c         fac1 = -6.d0*gmc2*mass(1)*ir3h(i)*irh(i)
c         axh(i) = fac1*xh(i)
c         ayh(i) = fac1*yh(i)
c         azh(i) = fac1*zh(i)  

c... formula of Venturini & Gallardo 2009
c         gm = mass(1) + mass(i)
c         call orbel_xv2aeq(xh(i),yh(i),zh(i),vxh(i),vyh(i),vzh(i),
c     &     gm,ialpha,a,e,q)
c         fac2 = a*dsqrt(1.d0-e*e)
c         fac1 = 3.d0*gmc2*mass(1)*irh(i)/fac2/fac2/fac2
c         axh(i) = fac1*xh(i)
c         ayh(i) = fac1*yh(i)
c         azh(i) = fac1*zh(i)  

         endif       
      enddo
!$OMP END PARALLEL DO

      return
      end      ! getacch_gr

c*************************************************************************
c                            DISCARD_MASSIVE_HYBRID.F
c*************************************************************************
c This subroutine checks to see if a terrestrial body should be discarded or
c merged.
c
c             Input:
c                 time          ==>  current time (real scalar)
c                 dt            ==>  time step  (real scalar)
c                 nbod          ==>  number of terrestrials bodies (int scalar)
c                 mass          ==>  mass of terrestrials (real array)
c                 xh,yh,zh      ==>  terrestrials position in helio coord 
c                                    (real arrays)
c                 vxh,vyh,vzh   ==>  terrestrials vel in helio coord 
c                                    (real arrays)
c                 njov          ==>  number of jovians bodies (int scalar)
c                 massj         ==>  mass of jovians (real array)
c                 xhj,yhj,zhj   ==> jovians position in helio coord 
c                                    (real arrays)
c                 vxhj,vyhj,vzhj ==> jovians vel in helio coord 
c                                    (real arrays)
c                 rmin,rmax      ==>  maximum and min distance from Sun
c                                     if <0  then don't check
c                                        (real scalar)
c                 rmaxu          ==>  maximum distance from Sun in not bound
c                                     if <0  then don't check
c                                        (real scalar)
c                 qmin           ==> Smallest perihelion distance
c                                      if <0  then don't check
c                                          (real scalar)
c                 lclose        ==> .true. --> marge terrestrials if they get too close. Read in that 
c                                    distance in io_init_pl
c                                          --> discard terrestrial if encounter with jovian
c                                      (logical*2 scalar)
c                 rpl           ==>  physical size of a planet.
c                                    (real array)
c                 rhill         ==>  size of a terrestrial hill's sphere.
c                                    (real array)
c                 isenc         ==>  0 --> No encounter during last dt
c                                    1 --> There was encounters
c                                     (integer scalar)
c                 eoff          ==> Amount of energy lost due to discards
c                                          (real scalar)
c                 mergelst      ==>  list of mergers (int array)
c                 mergecnt      ==>  count of mergers (int array)
c                 iecnt         ==>  Number of encounters (int*2 array)
c             Output:
c                 nbod          ==>  recalculated number of terrestrials bodies 
c                                       (int scalar)
c                 mass          ==>  recalculated mass of terrestrials bodies (real array)
c                 xh,yh,zh      ==>  recalculated position in helio coord 
c                                    (real arrays)
c                 vxh,vyh,vzh   ==>  recalculated vel in helio coord 
c                                    (real arrays)
c                 rpl           ==> recalculated physical sizes of terrestrials.
c                                    (real array)
c                 rhill         ==>  reordered size of terrestrials hill's sphere.
c                                    (real array)
c                 eoff          ==> Updated amount of energy lost from discards
c                                          (real scalar)
c
c
c Remarks: Based on discard_massive5. Adds check to see if a terrestrial has a close 
c          encounter with a jovian
c          iwhy =  0 no discard
c               =  n encounter with jovian n
c               = -1 planet to close to Sun
c               = -2 planet unbound from barycenter
c               = -3 planet too far from Sun
c               = -4 perihelion too small


      subroutine discard_massive_hybrid(time,dt,nbod,mass,xh,yh,zh,
     &     vxh,vyh,vzh,njov,massj,xhj,yhj,zhj,vxhj,vyhj,vzhj,
     &     rmin,rmax,rmaxu,qmin,lclose,rpl,
     &     rhill,isenc,mergelst,mergecnt,iecnt,eoff)

      include './swift.inc'

c...  Inputs: 
      real*8 time,dt
      integer isenc
      real*8 rmin,rmax,rmaxu,qmin
      logical*2 lclose
      integer mergelst(2,NTPMAX),mergecnt
      integer iecnt(NTPMAX)
      integer njov
      real*8 massj(njov),xhj(njov),yhj(njov),zhj(njov)
      real*8 vxhj(njov),vyhj(njov),vzhj(njov)


c...  Input and Output
      integer nbod
      real*8 mass(nbod),xh(nbod),yh(nbod),zh(nbod)
      real*8 vxh(nbod),vyh(nbod),vzh(nbod)
      real*8 eoff,rpl(nbod),rhill(nbod)

c...  internal
      integer iwhy(NTPMAX),i,iu,iflg,i1,i2,j,isperih(NTPMAX)
      real*8 xb(NTPMAX),yb(NTPMAX),zb(NTPMAX),rhj(NTPMAX)
      real*8 vxb(NTPMAX),vyb(NTPMAX),vzb(NTPMAX)
      real*8 rmin2,rmax2,rmaxu2,energy
      real*8 ei,ef,ke,pot,eltot(3),vdotr
      real*8 rh2,rb2,vb2,msys
      integer aflag
      real*8 r2,r2crit,r2crito,xr,yr,zr,vxr,vyr,vzr
      logical*1 lrflg(NTPMAX)
      character*1 cdummy

      save isperih

c-----
c...  Executable code 

c.... check for duplicate mergers
      do i=1,nbod
         lrflg(i) = .true.
      enddo
      do i=1,mergecnt
         i2 = mergelst(2,i)
         if(lrflg(i2)) then
            lrflg(i2) = .false.
         else
            mergelst(2,i) = -1
         endif
      enddo

c.... take care of mergers
      do i=1,mergecnt 
         i1 = mergelst(1,i)
         i2 = mergelst(2,i)
         vdotr = xh(i1)*vxh(i1)+yh(i1)*vyh(i1)+zh(i1)*vzh(i1)
         if (vdotr .gt. 0.d0) then
            isperih(i1) = 1
         else 
            isperih(i1) =-1
         endif
         if(i2.gt.0) then
            call discard_mass_reorder5(i2,nbod,mass,xh,yh,zh,
     &           vxh,vyh,vzh,rpl,rhill,isperih)
            do j=i+1,mergecnt
               if(mergelst(1,j).gt.i2) then
                  mergelst(1,j) = mergelst(1,j) - 1
               endif
               if(mergelst(2,j).gt.i2) then
                  mergelst(2,j) = mergelst(2,j) - 1
               endif
            enddo
         endif
      enddo

c...  set things up
      do i=1,nbod
         iwhy(i) = 0
      enddo

c...  check for position
      if( (rmin.ge.0.0) .or. (rmax.ge.0.0) .or. (rmaxu.ge.0.0) ) then
         rmin2 = rmin*rmin
         rmax2 = rmax*rmax
         rmaxu2 = rmaxu*rmaxu

         call coord_h2b_ter(nbod,mass,xh,yh,zh,vxh,vyh,vzh,
     &      njov,massj,xhj,yhj,zhj,vxhj,vyhj,vzhj,
     &      xb,yb,zb,vxb,vyb,vzb,msys)                 ! modify to properly account for the barycenter

!$OMP    PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) 
!$OMP&      PRIVATE(i,rh2,rb2,vb2,energy) 
!$OMP&      SHARED(nbod,xh,yh,zh,vxh,vyh,vzh,xb,yb,zb,vxb,vyb,vzb,msys,
!$OMP&            rmax,rmax2,rmin,rmin2,rmaxu,rmaxu2,iecnt,time,iwhy)
         do i=2,nbod
            rh2 = xh(i)**2 + yh(i)**2 + zh(i)**2
            if( (rmax.ge.0.0) .and. (rh2.gt.rmax2) ) then
!$OMP          CRITICAL
               write(*,*) rmax2,rh2,i
               write(*,*) 'Planet',i,' too far from Sun at t=',
     &              time
c               call flush()
               iwhy(i) = -3
!$OMP          END CRITICAL
            endif
            if( (rmin.ge.0.0) .and. (rh2.lt.rmin2) ) then
!$OMP          CRITICAL
               write(*,*) 'Planet',i,' too close from Sun at t=',
     &              time
c               call flush()
               iwhy(i) = -1
!$OMP          END CRITICAL
            endif
            if((iecnt(i).eq.0).and.(rmaxu.ge.0.0).and.
     &           (iwhy(i).eq.0)) then
               rb2 = xb(i)**2 + yb(i)**2 + zb(i)**2
               vb2 = vxb(i)**2 + vyb(i)**2 + vzb(i)**2
               energy = 0.5*vb2 - msys/sqrt(rb2)
               if( (energy.gt.0.0) .and. (rb2.gt.rmaxu2) ) then
!$OMP             CRITICAL
                  write(*,*) 'Planet',i,' is unbound and too far ',
     &                 'from barycenter at t=',time
c                  call flush()
                  iwhy(i) = -2
!$OMP             END CRITICAL
               endif
            endif
         enddo
!$OMP    END PARALLEL DO
      endif

c...  check perihelion distance
      if(qmin.ge.0.0) then
         call discard_mass_peri(time,nbod,iecnt,mass,xh,yh,zh,
     &       vxh,vyh,vzh,qmin,iwhy,isperih)
c         call flush()
      endif

c... check for terrestrial removal due to encounter with jovian
      if (lclose) then
         do j=2,3   !njov this is corrected to avoid ice 3
            call util_hills1(massj(1),massj(j),xhj(j),yhj(j),zhj(j),
     &        vxhj(j),vyhj(j),vzhj(j),rhj(j)) 
!$OMP       PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) 
!$OMP&         PRIVATE(i,xr,yr,zr,vxr,vyr,vzr,r2,r2crit,r2crito,aflag) 
!$OMP&         SHARED(j,nbod,xh,yh,zh,vxh,vyh,vzh,rhill,time,
!$OMP&                xhj,yhj,zhj,vxhj,vyhj,vzhj,rhj,dt,iwhy)
            do i=2,nbod
               if (iwhy(i).eq.0) then
                  xr = xh(i) - xhj(j)
                  yr = yh(i) - yhj(j)
                  zr = zh(i) - zhj(j)
                  vxr = vxh(i) - vxhj(j)
                  vyr = vyh(i) - vyhj(j)
                  vzr = vzh(i) - vzhj(j)
                  r2 = xr**2 + yr**2 + zr**2
                  r2crit = ( (rhill(i)+rhj(j))*1.d0 )**2
                  r2crito = r2crit*2.5d0
                  call rmvs_chk_ind(xr,yr,zr,vxr,vyr,vzr,dt,
     &                     r2crito,r2crit,aflag)
                  if(aflag.ne.0) then
!$OMP                CRITICAL
                     write(*,*) 'Planet',i,' too close to Jovian',
     &                   j,' at t=',time
                     write(*,*) sqrt(r2crit),sqrt(r2crito),sqrt(r2)
c                     call flush()
                     iwhy(i) = j+abs(aflag)-1
!$OMP                END CRITICAL
                  endif
               endif
            enddo
!$OMP       END PARALLEL DO
         enddo
      endif

c... discard things and writes out
      iu = 40
      i = 2
      iflg = 0
      do while(i.le.nbod) 
         if(iwhy(i).ne.0) then
            if(iflg.eq.0) then
               iflg = 1
               call anal_energy(nbod,mass,0.0d0,0.0d0,xh,yh,zh,
     &           vxh,vyh,vzh,ke,pot,ei,eltot)
            endif
            call io_discard_mass(1,time,i,mass(i),rpl(i),xh(i),yh(i),
     &           zh(i),vxh(i),vyh(i),vzh(i),iu,iwhy(i),cdummy)
            do j=i,nbod-1
               iwhy(j) = iwhy(j+1)
            enddo
            call discard_mass_reorder5(i,nbod,mass,xh,yh,zh,
     &           vxh,vyh,vzh,rpl,rhill,isperih)
         else
            i = i + 1
         endif
      enddo

      if(iflg.ne.0) then
         call anal_energy(nbod,mass,0.0d0,0.0d0,xh,yh,zh,
     &        vxh,vyh,vzh,ke,pot,ef,eltot)
         eoff = eoff + ei - ef
      endif

      return
      end       

c***********************************************************************
c                          COORD_H2B_TER.F
c***********************************************************************
c     PURPOSE: Converts from Heliocentric to Barycentric coords. 
c              Terrestrials only
c     ARGUMENTS:  Input is 
c                    nbod => number of terrestrials (integer)
c                    mass(*) =>  masses of terrestrials (real array)
c                    njov => number of jovians (integer)
c                    massj(*) =>  masses of jovians (real array)
c                    xh(*),yh(*),zh(*) => Terrestrials helio position (real array)
c                    xhj(*),yhj(*),zhj(*) => Jovians helio position (real array)
c                    vxh(*),vyh(*),vzh(*) => Terrestrials helio velocities (real array)
c                    vxhj(*),vyhj(*),vzhj(*) => Jovians helio velocities (real array)
c                 Returned are
c                    xb(*),yb(*),zb(*) => Terrestrial bary position (real array)
c                    vxb(*),vyb(*),vzb(*) => Terrestrial bary velocities (real array)
c                    msys => total mass of the system (real)
c
c Remarks: Adapted from coord_h2b

      subroutine coord_h2b_ter(nbod,mass,xh,yh,zh,vxh,vyh,vzh,
     &      njov,massj,xhj,yhj,zhj,vxhj,vyhj,vzhj,
     &      xb,yb,zb,vxb,vyb,vzb,msys)

      include './swift.inc'

c...  Inputs: 
      integer nbod,njov
      real*8 mass(nbod),massj(njov)
      real*8 xh(nbod),yh(nbod),zh(nbod)
      real*8 xhj(njov),yhj(njov),zhj(njov)
      real*8 vxh(nbod),vyh(nbod),vzh(nbod)
      real*8 vxhj(njov),vyhj(njov),vzhj(njov)

c...  Outputs:
      real*8 xb(nbod),yb(nbod),zb(nbod),msys
      real*8 vxb(nbod),vyb(nbod),vzb(nbod)

c...  Internals:
      integer n
      real*8 msysj
      real*8 xtmp,ytmp,ztmp,vxtmp,vytmp,vztmp

c----
c...  Executable code 

c...  First add up the jovians
      msysj = mass(1) + massj(2)
      xb(1) = massj(2)*xhj(2)
      yb(1) = massj(2)*yhj(2)
      zb(1) = massj(2)*zhj(2)
      vxb(1) = massj(2)*vxhj(2)
      vyb(1) = massj(2)*vyhj(2)
      vzb(1) = massj(2)*vzhj(2)
      do n=3,njov
         msysj = msysj + massj(n)
         xb(1) = xb(1) + massj(n)*xhj(n)
         yb(1) = yb(1) + massj(n)*yhj(n)
         zb(1) = zb(1) + massj(n)*zhj(n)
         vxb(1) = vxb(1) + massj(n)*vxhj(n)
         vyb(1) = vyb(1) + massj(n)*vyhj(n)
         vzb(1) = vzb(1) + massj(n)*vzhj(n)
      enddo

c...  Then add up the terrestrials
      msys = 0.0d0
      xtmp = 0.0d0
      ytmp = 0.0d0
      ztmp = 0.0d0
      vxtmp = 0.0d0
      vytmp = 0.0d0
      vztmp = 0.0d0
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) 
!$OMP&    PRIVATE(n) SHARED(nbod,mass,xh,yh,zh,vxh,vyh,vzh) 
!$OMP&    REDUCTION(+:msys,xtmp,ytmp,ztmp,vxtmp,vytmp,vztmp) 
      do n=2,nbod
         msys = msys + mass(n)
         xtmp = xtmp + mass(n)*xh(n)
         ytmp = ytmp + mass(n)*yh(n)
         ztmp = ztmp + mass(n)*zh(n)
         vxtmp = vxtmp + mass(n)*vxh(n)
         vytmp = vytmp + mass(n)*vyh(n)
         vztmp = vztmp + mass(n)*vzh(n)
      enddo
!$OMP END PARALLEL DO
      msys = msys + msysj
      xb(1) = xb(1) + xtmp 
      yb(1) = yb(1) + ytmp 
      zb(1) = zb(1) + ztmp 
      vxb(1) = vxb(1) + vxtmp 
      vyb(1) = vyb(1) + vytmp 
      vzb(1) = vzb(1) + vztmp 

c...  Compute the barycentric coords of the Sun
      xb(1) = -xb(1)/msys
      yb(1) = -yb(1)/msys
      zb(1) = -zb(1)/msys
      vxb(1) = -vxb(1)/msys
      vyb(1) = -vyb(1)/msys
      vzb(1) = -vzb(1)/msys

c...  Barycentric coords of terrestrials
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) PRIVATE(n) 
!$OMP&  SHARED(nbod,xb,yb,zb,vxb,vyb,vzb,xh,yh,zh,vxh,vyh,vzh) 
      do n=2,nbod
        xb(n) = xh(n) + xb(1)
        yb(n) = yh(n) + yb(1)
        zb(n) = zh(n) + zb(1)
        vxb(n) = vxh(n) + vxb(1)
        vyb(n) = vyh(n) + vyb(1)
        vzb(n) = vzh(n) + vzb(1)
      enddo
!$OMP END PARALLEL DO

      return
      end     ! coord_h2b_ter

c*************************************************************************
c                            IO_READ_LINE_R
c*************************************************************************
c read one line from real*4 binary file.
c
c      Input:
c            iu       ==> unit number to write to
c      Output:
C	     a        ==> semi-major axis or pericentric distance if a parabola
c                          (real scalar)
c            e        ==> eccentricity (real scalar)
C            inc      ==> inclination  (real scalar)
C            capom    ==> longitude of ascending node (real scalar)
C	     omega    ==> argument of perihelion (real scalar)
C	     capm     ==> mean anomoly(real scalar)
c       Returns:
c      io_read_line_r    ==>   =0 read ok
c                           !=0 read failed is set to iostat variable
c
c Remarks: 
c Authors:  Hal Levison 
c Date:    2/22/94
c Last revision: 

      integer function io_read_line_r(iu,id,a,e,inc,capom,omega,capm) 

      include './swift.inc'
      include './io.inc'

c...  Inputs: 
      integer iu

c...  Output: 
      integer id
      real*8 a,e,inc,capom,omega,capm

c...  Internals
      integer*2 id2
      real*4 a4,e4,inc4,capom4,omega4,capm4
      integer ierr

c----
c...  Executable code 

      read(iu,iostat=ierr) id2,a4,e4,inc4,capom4,omega4,capm4
      io_read_line_r = ierr
      if(ierr.ne.0) then
         return
      endif

      id = id2

      a = a4
      e = e4
      inc = inc4
      capom = capom4
      capm = capm4
      omega = omega4

      return
      end      ! io_read_line_r
c--------------------------------------------------------------------------

c*************************************************************************
c                            IO_READ_HDR_R
c*************************************************************************
c read in header part of the real*4 file
c
c             Input:
c                 iu            ==> unit number to write to
c             Output:
c                 time          ==>  current time (real scalar)
c                 nbod          ==>  number of massive bodies (int scalar)
c                 nleft         ==>  number of active tp (int scalar)
c
c             Returns:
c               io_read_hdr_r     ==>   =0 read ok
c                                    !=0 read failed is set to iostat variable
c Remarks: 
c Authors:  Hal Levison 
c Date:    2/22/94
c Last revision: 

      integer function io_read_hdr_r(iu,time,nbod,nleft) 

      include './swift.inc'
      include './io.inc'

c...  Inputs: 
      integer iu

c...  Output
      integer nbod,nleft
      real*8 time

c...  Internals
      real*8 ttmp
c      real*4 ttmp
      integer*2 nleft2,nbod2
      integer ierr

c----
c...  Executable code 

      read(iu,iostat=ierr) ttmp,nbod2,nleft2
      io_read_hdr_r = ierr
      if(ierr.ne.0) then
         return
      endif

      nbod = nbod2
      nleft = nleft2
      time = ttmp

      return
      end     ! io_read_hdr_r.f

c*************************************************************************
c                            IO_READ_MASS_R
c*************************************************************************
c read in the mass file.
c
c             Output:
c                 time          ==>  current time (real scalar)
c                 nbod          ==>  number of massive bodies (int scalar)
c                 mass          ==>  mass of bodies (real array)
c                 iu              ==> unit number to read to
c
c             Returns:
c               io_read_mass     ==>   =0 read ok
c                                    !=0 read failed is set to iostat variable
c
c Remarks: Based on io_read_frame
c Authors:  Hal Levison 
c Date:    1/9/97
c Last revision: 11/2/99

      integer function io_read_mass_r(time,nbod,mass,iu)

      include './swift.inc'
      include './io.inc'

c...  Inputs: 
      integer iu

c...  Outputs
      integer nbod
      real*8 mass(nbod),time

c...  Internals
      real*4 mass4(NTPMAX)
      real*8 ttmp
c      real*4 ttmp
      integer*2 nbod2
      integer i,ierr

c----
c...  Executable code 

      read(iu,iostat=ierr) ttmp,nbod2
      io_read_mass_r = ierr
      if(ierr.ne.0) then
         return
      endif

      read(iu,iostat=ierr) (mass4(i),i=1,nbod2)
      io_read_mass_r = ierr
      if(ierr.ne.0) then
         return
      endif

      nbod=nbod2
      do i=1,nbod
         mass(i) = mass4(i)
      enddo
      time = ttmp

      return
      end      ! io_read_mass_r

c*************************************************************************
c                            IO_DISCARD_MASS
c*************************************************************************
c Write out information about a discarded massive body.
c
c             Input:
c                 init          ==>  initiize flag if = 0 initialize and return
c                                                     = 1 run through 
c                 id            ==> particle number (int scalar)
c                 time          ==>  current time (real scalar)
c                 m1            ==>  Mass of pl (real scalar)
c                 r1            ==>  Radius of pl 2 (real scalar)
c                 x1,y1,z1      ==>  current position of pl 1 in helio coord 
c                                    (real scalar)
c                 vx1,vy1,vz1   ==>  current velocity of pl 1 in helio coord 
c                                    (real scalar)
c                 iu            ==> IO unit (int scalar)
c                 iwhy          ==> reason for discard (int scalar)
c                 fopenstat     ==>  The status flag for the open 
c                                      statements of the output files.  
c                                          (character*80)
c Remarks: 
c Authors:  Hal Levison 
c Date:    12/30/96
c Last revision: 9/11/09

      subroutine io_discard_mass(init,time,id,m1,r1,x1,y1,z1,vx1,vy1,
     &     vz1,iu,iwhy,fopenstat)

      include './swift.inc'
      include './io.inc'

c...  Inputs: 
      integer iwhy,iu,init,id
      real*8 time
      real*8 m1,r1
      real*8 x1,y1,z1
      real*8 vx1,vy1,vz1
      character*(*) fopenstat

      character*80 filediscard
      common /discard/ filediscard

c...  Internals
      integer ierr

c----
c...  Executable code 

      if(init.eq.0) then

         call io_open(iu,filediscard,fopenstat,
     &        'FORMATTED',ierr)

c...     if there was an error and fopenstat='append' then
c...     try to open as new
         if(ierr.ne.0) then  
            if( (fopenstat(1:6).eq.'append') .or. 
     &           (fopenstat(1:6).eq.'APPEND') ) then
               call io_open(iu,filediscard,'new','FORMATTED',
     &              ierr)
            endif
         endif

         if(ierr.ne.0) then
            write(*,*) ' SWIFT ERROR: in io_discard_mass: '
            write(*,*) '    Could not open discard output file'
            call util_exit(1)
         endif
         close(unit = iu)
         return   ! <=== NOTE!!!!
      else
         call io_open(iu,filediscard,'append','FORMATTED',
     &        ierr)
      endif

      write(iu,1000) time,iwhy
 1000 format(1x,1p1e23.16,1x,i4)

      write(iu,2000) id,m1,r1
 2000 format('-1',1x,i5,1x,2(1p1e23.16,1x))
      write(iu,3000) x1,y1,z1
 3000 format(3(1p1e23.16,1x))
      write(iu,3000) vx1,vy1,vz1

      close(unit = iu)
      return
      end                       ! io_discard_mass.f
c--------------------------------------------------------------------------

c*************************************************************************
c                            IO_DISCARD_MERGE
c*************************************************************************
c Write out information about a merger.
c
c             Input:
c                 time          ==>  current time (real scalar)
c                 ip1,ip2       ==>  planets to merge (real scalar)
c                 m1            ==>  Mass of pl 1 (real scalar)
c                 r1            ==>  Radius of pl 1 (real scalar)
c                 x1,y1,z1      ==>  current position of pl 1 in helio coord 
c                                    (real arrays)
c                 vx1,vy1,vz1   ==>  current velocity of pl 1 in helio coord 
c                                    (real arrays)
c                 m2            ==>  Mass of pl 2 (real scalar)
c                 r2            ==>  Radius of pl 2 (real scalar)
c                 x2,y2,z2      ==>  current position of pl 2 in helio coord 
c                                    (real arrays)
c                 vx2,vy2,vz2   ==>  current velocity of pl 2 in helio coord 
c                                    (real arrays)
c                 mn            ==>  Mass of new pl  (real scalar)
c                 rn            ==>  Radius of new pl (real scalar)
c                 xn,yn,zn      ==>  current position of new pl in helio coord 
c                                    (real arrays)
c                 vxn,vyn,vzn   ==>  current velocity of new pl in helio coord 
c                                    (real arrays)
c                 nleft           ==>  number of active test bodies(int scalar)
c
c Remarks: 
c Authors:  Hal Levison 
c Date:    12/30/96
c Last revision: 

      subroutine io_discard_merge(time,ip1,ip2,m1,r1,x1,y1,z1,vx1,vy1,
     &     vz1,m2,r2,x2,y2,z2,vx2,vy2,vz2,mn,rn,xn,yn,zn,vxn,vyn,vzn)

      include './swift.inc'
      include './io.inc'

c...  Inputs: 
      integer ip1,ip2
      real*8 time
      real*8 m1,r1
      real*8 x1,y1,z1
      real*8 vx1,vy1,vz1
      real*8 m2,r2
      real*8 x2,y2,z2
      real*8 vx2,vy2,vz2
      real*8 mn,rn
      real*8 xn,yn,zn
      real*8 vxn,vyn,vzn

      character*80 filediscard
      common /discard/ filediscard
      

c...  Internals
      integer ierr,iu

c----
c...  Executable code 

      iu = 40

      call io_open(iu,filediscard,'append','FORMATTED',ierr)

      write(iu,1000) time
 1000 format(1x,1p1e23.16,'  2')

      write(iu,2000) ip1,m1,r1
 2000 format('-1',1x,i5,1x,2(1p1e23.16,1x))
      write(iu,3000) x1,y1,z1
 3000 format(3(1p1e23.16,1x))
      write(iu,3000) vx1,vy1,vz1

      write(iu,2000) ip2,m2,r2
      write(iu,3000) x2,y2,z2
      write(iu,3000) vx2,vy2,vz2

      write(iu,4000) ip1,mn,rn
 4000 format('+1',1x,i5,1x,2(1p1e23.16,1x))
      write(iu,3000) xn,yn,zn
      write(iu,3000) vxn,vyn,vzn

      close(unit = iu)
      return
      end                       ! io_discard_merge.f
c--------------------------------------------------------------------------
